In this current project, we have 3 groups of mice. A control group that was fed 10% fat diet, an oleate group which was fed with an oleate diet and a HFD group which was fed with 60% fat diet. The project included 19 mice. We collected feces, contents from the ileum, cecum and colon and we performed 16S rRNA sequencing on them. The DNA extraction and sequencing was done in the university of Michigan Microbiome core. The details of sequencing procedures are present in the Microbiome_core_sequencing_details.xlsx excel file.
First, we created a file that assign the fastq files to the samples. The file is called oleate.files in the data directory.
We used mothur to preprocess the data and to clasify OTUs. An executable version of mothur for windows (version 1.44.3) was installed in the folder containing the fastq files. The original fastq files were processed using mothur pipeline using the following code:
# setting the input and output directories
set.dir(input=Data, output= ButyrateResults)
# making contigs
make.contigs(file=oleate.files, processors=4)
screen.seqs(fasta=current, group=current, maxambig=0, maxlength=275)
unique.seqs(fasta=current)
count.seqs(name=current, group=current)
summary.seqs(count=current)
# make aan align reference file from silva original reference file. The silva reference files (version 138.1) were downloaded from <https://mothur.org/wiki/silva_reference_files/>
pcr.seqs(fasta=silva.nr_v138_1.align, start=11894, end=25319, keepdots=F, processors=8)
# align sequences
align.seqs(fasta= oleate.trim.contigs.good.unique.fasta, reference= silva.nr_v138_1.pcr.align)
summary.seqs(fasta= current, count= current)
screen.seqs(fasta=current, count=current, summary=current, start=1968, end=11550, maxhomop=8)
summary.seqs(fasta= current, count= current)
filter.seqs(fasta=current, vertical=T, trump=.)
unique.seqs(fasta=current, count=current)
pre.cluster(fasta=current, count=current, diffs=2)
chimera.vsearch(fasta=current, count=current, dereplicate=t)
remove.seqs(fasta=current, accnos=current)
summary.seqs(fasta=current, count=current)
# classify sequences
classify.seqs(fasta=current, count=current, reference= silva.nr_v138_1.pcr.align, taxonomy= silva.nr_v138_1.tax, cutoff=80)
remove.lineage(fasta=current, count=current, taxonomy=current, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota)
summary.tax(taxonomy=current, count=current)
rename.file(taxonomy=oleate.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.0.03.cons.taxonomy, shared= oleate.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.shared)
# investigate the read counts in each sample and selecting the least count for subsampling
count.groups(shared=oleate.opti_mcc.shared)
sub.sample(shared= current, size= 2766)
rarefaction.single(shared=current, calc=sobs, freq=100)
summary.single(shared=current, calc=nseqs-coverage-sobs-invsimpson, subsample=T)
# calculating distances and principal coordinate analysis
dist.shared(shared=current, calc=thetayc-jclass, subsample=t)
pcoa(phylip=current)
nmds(phylip=current)
nmds(phylip=current, mindim=3, maxdim=3)
# we created a set of txt files in excel to reflect the design and metadata of the samples. these files are present in the data directory under the names `condition.design`, `site.design` and `oleate.metadata`. Next, we proceeded with performing AMOVA to detect significance
amova(phylip=current, design=condition.design)
amova(phylip=current, design=site.design)
homova(phylip=current, design=condition.design)
homova(phylip=current, design=site.design)
We created an R project and copied the following files to the project directory:
1- oleate.opti_mcc.shared (the OTU table of all samples)
2- oleate.opti_mcc.0.03.subsample.shared (the OTU table after subsampling)
3-oleate.metadata (a txt file that includes the sample metadata)
4- oleate.taxonomy (the taxonomy file that was generated by mothur. we modified this file in excel. As it is from mothur, there are not column headers for order, family genus etc. in the cons.taxonomy file and this must be fixed first. This is how I did it: read cons.taxonomy file into excel and choose semicolon as a separator so each tax level is therefore in its own column. then delete header ‘taxonomy’ and put in appropriate header for each column (order or family or genus etc. Then select all, press ctrl+shift+H on mac to find and replace. Replate all (*) with nothing)
## Install and load required libraries
library("phyloseq")
library("ggplot2")
library("plyr")
library("vegan")
library("grid")
library("directlabels")
library("knitr")
library("clustsig")
library("ellipse")
library("ape")
library("DESeq2")
library("microbial")
library("VennDiagram")
library("microbiome")
library("microbiomeutilities")
## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'
library("viridis")
library("RColorBrewer")
library("ggpubr")
library("patchwork")
## import data
sharedFile = read.table("oleate.opti_mcc.shared")
sharedFile = t(sharedFile) ## transform data
rownames(sharedFile) = sharedFile[,1] ## define rowname
colnames(sharedFile) = sharedFile[2,] ## define column names
sharedFile <-as.matrix(sharedFile)
sharedFile = sharedFile[,2:57] ## only include columns that reflects OTU numbers (remember order [row,colum]) 2 -> 57 for me
## as I had 57 samples
sharedFile = sharedFile[4:819,] ## and which rows you want. I have 819 OTUs so I changed this number to 819
# original => sharedFile = sharedFile[4:819,]
class(sharedFile) <- "numeric"
head(sharedFile) ## look at head first 7 lines to see it's made correclty
## ce4406 ce4407 ce4408 ce4410 ce4411 ce4412 ce4413 ce4414 ce4415 ce4416
## Otu001 648 490 686 411 726 838 761 175 144 1520
## Otu002 1 143 1693 1037 794 548 356 418 149 1075
## Otu003 867 743 1670 733 1171 559 953 165 202 564
## Otu004 310 8 1379 808 903 419 1487 953 46 192
## Otu005 483 704 1357 2703 825 735 381 416 188 2496
## Otu006 0 0 0 1 170 2 4 49 40 365
## ce4417 ce4418 ce4419 ce4420 ce4421 ce4422 ce4423 ce4424 co4406 co4407
## Otu001 1118 1091 1691 703 936 1011 54 1452 1170 257
## Otu002 2412 6002 1290 586 1216 1117 2346 48 5 843
## Otu003 1201 1204 871 980 890 922 1522 388 406 1539
## Otu004 29 291 1469 21 6 4 3 20 1066 14
## Otu005 2134 969 1805 1737 1385 1602 1220 391 292 188
## Otu006 123 1317 97 245 237 272 146 2061 3 2
## co4408 co4410 co4411 co4412 co4413 co4414 co4415 co4416 co4417 co4418
## Otu001 1823 856 397 705 356 392 45 201 4331 356
## Otu002 767 177 220 404 24 442 261 1187 1901 1269
## Otu003 1116 517 909 996 394 499 490 498 664 1311
## Otu004 1089 564 518 468 660 1151 27 133 481 13
## Otu005 264 400 795 778 310 875 168 31 804 483
## Otu006 0 0 59 1 1 110 19 195 1149 66
## co4419 co4420 co4421 co4422 co4423 co4424 co4428 f4406 f4407 f4408 f4410
## Otu001 182 188 1372 944 1621 233 1166 7021 9127 2210 3181
## Otu002 2380 1814 815 390 778 2231 10 31 1790 14 173
## Otu003 567 499 527 791 1318 1249 225 1108 1447 679 726
## Otu004 31 484 35 0 18 4 6 1221 2604 386 2143
## Otu005 281 212 501 519 584 237 148 1300 380 1279 628
## Otu006 240 39 280 237 667 376 1572 3 0 0 1
## f4411 f4412 f4413 f4414 f4415 f4416 f4417 f4418 f4419 f4420 f4421 f4422
## Otu001 6323 8943 2011 4584 4530 7401 3309 6001 7938 4862 1619 1334
## Otu002 1367 30 99 2063 1193 898 472 265 3751 3667 1244 1022
## Otu003 1242 1133 1532 1163 687 902 536 443 1878 554 1781 1668
## Otu004 2441 716 1082 4784 4027 3709 3814 2177 1569 2593 80 21
## Otu005 1154 1847 1024 816 591 1263 1511 521 1505 1013 272 558
## Otu006 786 1 0 1648 862 1234 1175 931 2214 1074 2252 1979
## f4423 f4424 f4428
## Otu001 3749 11842 791
## Otu002 1421 1194 996
## Otu003 610 2459 971
## Otu004 147 23 5
## Otu005 1367 64 324
## Otu006 798 506 1765
## Import subsampled otu matrix (26880 seqs)
sharedsubFile = read.table('oleate.opti_mcc.0.03.subsample.shared')
sharedsubFile = t(sharedsubFile)
rownames(sharedsubFile) = sharedsubFile[,1]
colnames(sharedsubFile) = sharedsubFile[2,]
sharedsubFile = sharedsubFile[,2:57]
sharedsubFile = sharedsubFile[4:494,]
class(sharedsubFile) <- "numeric"
head(sharedsubFile)
## ce4406 ce4407 ce4408 ce4410 ce4411 ce4412 ce4413 ce4414 ce4415 ce4416
## Otu001 335 271 174 72 204 358 233 131 144 312
## Otu002 1 68 361 199 235 236 103 311 149 237
## Otu003 468 380 379 139 359 254 289 125 202 108
## Otu004 145 2 312 128 261 186 453 680 46 37
## Otu005 245 348 299 492 254 319 124 308 188 509
## Otu006 0 0 0 0 38 1 1 31 40 84
## ce4417 ce4418 ce4419 ce4420 ce4421 ce4422 ce4423 ce4424 co4406 co4407
## Otu001 187 161 311 109 212 197 10 137 456 145
## Otu002 414 825 219 86 255 194 590 3 2 528
## Otu003 213 166 156 157 197 176 380 30 160 916
## Otu004 6 38 284 3 0 1 0 1 407 9
## Otu005 412 134 311 262 300 295 299 34 128 109
## Otu006 20 173 20 46 47 44 30 184 2 2
## co4408 co4410 co4411 co4412 co4413 co4414 co4415 co4416 co4417 co4418
## Otu001 473 406 172 289 160 168 41 127 611 138
## Otu002 224 82 101 156 11 181 240 770 281 465
## Otu003 304 272 445 419 180 214 459 308 85 543
## Otu004 317 258 253 167 321 469 27 86 69 7
## Otu005 83 202 400 294 147 373 155 19 105 172
## Otu006 0 0 24 0 0 44 18 121 169 18
## co4419 co4420 co4421 co4422 co4423 co4424 co4428 f4406 f4407 f4408 f4410
## Otu001 84 80 358 335 340 72 237 933 1161 548 615
## Otu002 873 792 208 118 166 744 1 4 230 2 29
## Otu003 212 202 130 290 274 413 48 141 199 168 115
## Otu004 18 223 9 0 4 1 2 151 360 105 410
## Otu005 111 73 132 184 122 82 24 154 52 311 113
## Otu006 95 14 68 96 126 127 293 0 0 0 0
## f4411 f4412 f4413 f4414 f4415 f4416 f4417 f4418 f4419 f4420 f4421 f4422
## Otu001 929 1220 419 585 727 779 584 1196 788 533 250 228
## Otu002 196 5 22 259 167 105 74 49 390 416 224 144
## Otu003 171 167 296 144 124 95 99 87 195 75 271 261
## Otu004 384 99 204 639 655 415 673 421 145 281 15 5
## Otu005 179 283 209 99 101 125 256 93 145 131 40 87
## Otu006 106 0 0 177 135 128 212 197 229 116 332 338
## f4423 f4424 f4428
## Otu001 401 1571 135
## Otu002 158 159 161
## Otu003 65 348 150
## Otu004 12 3 0
## Otu005 137 11 59
## Otu006 95 66 275
dim(sharedsubFile)
## [1] 491 56
sharedsubFile <-as.matrix(sharedsubFile)
## Import taxonomy file
## Copy and paste into notepad and save. Then carry on:
taxFile = read.table('oleate.taxonomy', header=T, sep='\t')
rownames(taxFile) = taxFile[,1]
taxFile = taxFile[,3:8]
taxFile = as.matrix(taxFile)
head(taxFile)
## Kingdom Phylum Class Order
## Otu001 "Bacteria" "Firmicutes" "Bacilli" "Lactobacillales"
## Otu002 "Bacteria" "Verrucomicrobiota" "Verrucomicrobiae" "Verrucomicrobiales"
## Otu003 "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales"
## Otu004 "Bacteria" "Firmicutes" "Bacilli" "Lactobacillales"
## Otu005 "Bacteria" "Firmicutes" "Clostridia" "Lachnospirales"
## Otu006 "Bacteria" "Firmicutes" "Bacilli" "Erysipelotrichales"
## Family Genus
## Otu001 "Streptococcaceae" "Lactococcus"
## Otu002 "Akkermansiaceae" "Akkermansia"
## Otu003 "Bacteroidaceae" "Bacteroides"
## Otu004 "Lactobacillaceae" "Lactobacillus"
## Otu005 "Lachnospiraceae" "Lachnospiraceae"
## Otu006 "Erysipelotrichaceae" "Dubosiella"
## import metadata file
metaFile = read.table('oleate.metadata', header=T, sep='\t')
rownames(metaFile) = metaFile[,1]
metaFile = metaFile[,1:4]
head(metaFile)
## group site mouse condition
## ce4406 ce4406 cecum 4406 10%Fat
## ce4407 ce4407 cecum 4407 10%Fat
## ce4408 ce4408 cecum 4408 10%Fat
## ce4410 ce4410 cecum 4410 10%Fat
## ce4411 ce4411 cecum 4411 10%Fat
## ce4412 ce4412 cecum 4412 10%Fat
## Create phyloseq object
OTU = otu_table(sharedFile, taxa_are_rows = TRUE)
OTUsub = otu_table(sharedsubFile, taxa_are_rows = TRUE)
TAX = tax_table(taxFile)
META = sample_data(metaFile)
physeq = phyloseq(OTU, TAX, META)
physeqSub = phyloseq(OTUsub, TAX, META)
physeqSub
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 491 taxa and 56 samples ]
## sample_data() Sample Data: [ 56 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 491 taxa by 6 taxonomic ranks ]
## Get rid of any OTUs not present in any samples and get relative abundance
microSub <- prune_taxa(taxa_sums(physeqSub) > 0, physeqSub)
microSubRel = transform_sample_counts(microSub, function(x) x / sum(x) )
microSubRelFilt = filter_taxa(microSubRel, function(x) mean(x) > 1e-5, TRUE)
# for subsampled shared file
# create a tree and a new phyloseq object
random_tree = rtree(ntaxa(microSub), rooted=TRUE, tip.label=taxa_names(microSub))
plot(random_tree)
microSubt = merge_phyloseq(microSub, random_tree)
# Plotting rarefaction
rarecurve(t(otu_table(microSubt)), step=50, cex=0.5)
# Plot a prettier rarefaction curve
# set seed
set.seed(1)
subsamples <- seq(0, 5000, by=100)[-1]
#subsamples = c(10, 5000, 10000, 20000, 30000)
p <- plot_alpha_rcurve(microSub, index="observed",
subsamples = subsamples,
lower.conf = 0.025,
upper.conf = 0.975,
group="condition",
label.color = "brown3",
label.size = 3,
label.min = TRUE)
# change color of line
mycols <- c("brown3", "steelblue","grey50")
p <- p + scale_color_manual(values = mycols) +
scale_fill_manual(values = mycols)
print(p)
# Plotting taxonomy
plot_bar(microSubt, fill="Phylum") +
facet_wrap(~condition, scales = "free_x", nrow = 1) +
theme(text = element_text(size=16))+geom_bar(stat="identity") +
scale_fill_brewer(type = "seq", palette = "Spectral") +
theme(
legend.background = element_rect(
fill = "lemonchiffon",
colour = "grey50",
size = 1
)
)
plot_bar(microSubt, fill="Phylum") +
facet_wrap(~site, scales = "free_x", nrow = 1) +
theme(text = element_text(size=16))+geom_bar(stat="identity") +
scale_fill_brewer(type = "seq", palette = "Spectral") +
theme(
legend.background = element_rect(
fill = "lemonchiffon",
colour = "grey50",
size = 1
)
)
plot_bar(microSub, fill="Phylum") +
facet_grid(condition~site, scales = "free") + ## separate facets
theme(text = element_text(size=16)) + ## increase font of all elements
geom_bar(stat="identity")+
scale_fill_brewer(type = "seq", palette = "Spectral") +
theme(
legend.background = element_rect(
fill = "lemonchiffon",
colour = "grey50",
size = 1
)
)
# Plotting PCoA
my.PCoA2 <- ordinate(microSub, "PCoA", "bray")
plot_ordination(microSub, my.PCoA2, type = "group", color = "condition", shape = "site")+theme(text = element_text(size=20))+ geom_point(size=5) + scale_color_manual(values=c("brown3", "steelblue","grey50"))
plot_ordination(microSub, my.PCoA2, type = "group", color = "condition")+ facet_wrap(~site, scales = "free_x", nrow = 1) + geom_point(size=5)+theme(text = element_text(size=20)) + scale_color_manual(values=c("brown3", "steelblue","grey50"))
plot_ordination(microSub, my.PCoA2, type = "group", color = "condition", shape = "site")+theme(text = element_text(size=20))+ geom_point(size=5)+
geom_line() + scale_color_manual(values=c("brown3", "steelblue","grey50"))
plot_ordination(microSub, my.PCoA2, type = "group", color = "condition")+theme(text = element_text(size=20))+ geom_point(size=5)+ stat_ellipse(level=0.9) + scale_color_manual(values=c("brown3", "steelblue","grey50"))
# PCoA plot using the unweighted UniFrac as distance
wunifrac_dist = phyloseq::distance(microSubt, method="unifrac", weighted=F)
ordination = ordinate(microSubt, method="PCoA", distance=wunifrac_dist)
plot_ordination(microSubt, ordination, color="condition") + theme(aspect.ratio=1)+
ggtitle("PCoA: unweigthed Unifrac")+geom_point(size=5) + scale_color_manual(values=c("brown3", "steelblue","grey50"))
#Plot PCoA using the weighted UniFrac as distance
wunifrac_dist = phyloseq::distance(microSubt, method="unifrac", weighted=T)
ordUF = ordinate(microSubt, method="PCoA", distance=wunifrac_dist)
plot_ordination(microSubt, ordUF, color = "condition") +
ggtitle("PCoA: weigthed Unifrac")+geom_point(size=5) + scale_color_manual(values=c("brown3", "steelblue","grey50"))
## Normalization and plotting figures by group
# Normalization
phy <- normalize(microSubt, method = "relative")
## Normalization using relative method
# Plot taxonomy by group
plotbar(phy,level="Phylum", group="condition") + theme(axis.text.x = element_text(size = 16))
plotbar(phy,level="Phylum", group="site") + theme(axis.text.x = element_text(size = 16))
# Plot alpha diversity
plotalpha(microSubt, group = "condition", color = c("brown3","grey50", "steelblue"))
plotalpha(microSubt, group = "site")
# Test for significance between groups
beta_condition <-betatest(phy,group="condition")
## Do PERMANOVA for: condition
beta_site <-betatest(phy,group="site")
## Do PERMANOVA for: site
# Biomarkers discovery and plotting
bio <- biomarker(microSubt,group="condition")
## Normalization using relative method
plotmarker(bio,level="Genus")
# LefSe testing and plotting
lda <- ldamarker(microSubt,group="condition")
## Normalization using relative method
write.csv(lda, "lda.csv", row.names = FALSE)
# open the csv file in excel, separate the taxonomy to kingdom, phylum,... and save it as lda2.csv so that the names can appear better on the graph
lda2 <- read.csv(file = 'lda2.csv')
plotLDA(lda2,group=c("10%Fat","60%Fat"),lda=5,pvalue=0.05)+theme(axis.text.y = element_text(size = 16))+theme(text = element_text(size=20))
plotLDA(lda2,group=c("10%Fat","Oleate"),lda=5,pvalue=0.05)+theme(axis.text.y = element_text(size = 16))+theme(text = element_text(size=20))
plotLDA(lda2,group=c("Oleate","60%Fat"),lda=5,pvalue=0.05) +theme(axis.text.y = element_text(size = 16))+theme(text = element_text(size=20))
## Modified difftest function from microbial package to test for significance
# create difftest2 function
difftest2<-function(physeq,group,pvalue=0.05,padj=NULL,log2FC=0,gm_mean=TRUE,fitType="local",quiet=FALSE){
if(!taxa_are_rows(physeq)){
physeq<-t(physeq)
}
otu <- as(otu_table(physeq),"matrix")
tax <- as.data.frame(as.matrix(tax_table(physeq)))
colData<-as(sample_data(physeq),"data.frame")
colData$condition<-colData[,group]
contrasts<-levels(factor(unique(colData$condition)))[1:2]
if(isTRUE(gm_mean)){
countData<-round(otu, digits = 0)
}else{
countData<-round(otu, digits = 0)+1
}
dds <- DESeqDataSetFromMatrix(countData, colData, as.formula(~condition))
if(isTRUE(gm_mean)){
geoMeans = apply(counts(dds), 1, gm_mean)
dds = estimateSizeFactors(dds, geoMeans = geoMeans)
}
dds <- DESeq(dds, fitType=fitType)
res <- results(dds,contrast=c("condition",contrasts),cooksCutoff = FALSE)
res_tax = cbind(as.data.frame(res), as.matrix(countData[rownames(res), ]))
if(!is.null(padj)){
pval<-padj
sig <- rownames(subset(res,padj<pval&abs(log2FoldChange)>log2FC))
}else{
pval<-pvalue
sig <- rownames(subset(res,pvalue<pval&abs(log2FoldChange)>log2FC))
}
res_tax$Significant<- "No"
res_tax$Significant <- ifelse(rownames(res_tax) %in% sig, "Yes", "No")
res_tax <- cbind(res_tax, tax[rownames(res),])
return(as.data.frame(res_tax))
}
# Test for significant bacteria using microbial package
mic_res <- difftest2(microSubt,group="condition", gm_mean=FALSE)
## converting counts to integer mode
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## estimating size factors
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## final dispersion estimates
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## fitting model and testing
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## fitting model and testing
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
plotdiff(mic_res,level="Genus",padj=0.001,log2FC = 3,fontsize.y = 14)
## Test for significant OTUs using DESeq2
sample_data(microSubt)$condition <- as.factor(sample_data(microSubt)$condition)
ds = phyloseq_to_deseq2(microSubt, ~ condition)
## converting counts to integer mode
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
ds = DESeq(ds)
## estimating size factors
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## final dispersion estimates
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## fitting model and testing
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## fitting model and testing
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
alpha = 0.01
# compare Oleate and 60%Fat
res = results(ds, contrast=c("condition", "Oleate", "60%Fat"), alpha=alpha)
res = res[order(res$padj, na.last=NA), ]
res_sig = res[(res$padj < alpha), ]
res_sig
## log2 fold change (MLE): condition Oleate vs 60%Fat
## Wald test p-value: condition Oleate vs 60%Fat
## DataFrame with 25 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## Otu004 200.56502 6.13780 0.634629 9.67149 3.98541e-22 7.81140e-20
## Otu031 17.90550 3.11457 0.399384 7.79842 6.26851e-15 6.14314e-13
## Otu028 18.25563 4.51975 0.610052 7.40880 1.27449e-13 8.32670e-12
## Otu029 21.61460 -5.03417 0.761618 -6.60984 3.84734e-11 1.80406e-09
## Otu054 7.35816 2.95435 0.448766 6.58327 4.60220e-11 1.80406e-09
## ... ... ... ... ... ... ...
## Otu050 8.71757 -2.90784 0.699150 -4.15910 3.19497e-05 0.000298197
## Otu027 21.35176 2.17939 0.534261 4.07927 4.51775e-05 0.000402491
## Otu040 11.60680 1.84137 0.487587 3.77648 1.59059e-04 0.001355458
## Otu078 2.00088 -3.40305 0.941310 -3.61523 3.00082e-04 0.002450666
## Otu062 5.35700 1.88742 0.584235 3.23058 1.23540e-03 0.009685536
res_sig = cbind(as(res_sig, "data.frame"), as(tax_table(microSub)[rownames(res_sig), ], "matrix"))
ggplot(res_sig, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
geom_jitter(size=5, width = 0.2) + theme(text = element_text(size=16)) +
theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
# compare 10%Fat and 60%Fat
res2 = results(ds, contrast=c("condition", "10%Fat", "60%Fat"), alpha=alpha)
res2 = res2[order(res2$padj, na.last=NA), ]
res_sig2 = res2[(res2$padj < alpha), ]
res_sig2
## log2 fold change (MLE): condition 10%Fat vs 60%Fat
## Wald test p-value: condition 10%Fat vs 60%Fat
## DataFrame with 53 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## Otu007 92.6627 -8.96186 0.564355 -15.87981 8.74359e-57 1.30280e-54
## Otu008 60.6831 -5.59715 0.559648 -10.00120 1.50568e-23 7.47822e-22
## Otu036 13.0374 -7.57488 0.756046 -10.01907 1.25680e-23 7.47822e-22
## Otu033 12.2718 -7.32904 0.771735 -9.49683 2.16382e-21 8.06022e-20
## Otu029 21.6146 -7.82794 0.852785 -9.17927 4.34030e-20 1.12102e-18
## ... ... ... ... ... ... ...
## Otu020 28.00857 1.19042 0.404478 2.94311 0.00324931 0.00949308
## Otu018 30.12132 1.71064 0.576419 2.96771 0.00300032 0.00949308
## Otu050 8.71757 -2.02837 0.687046 -2.95231 0.00315407 0.00949308
## Otu013 64.09575 1.45752 0.497800 2.92791 0.00341245 0.00977797
## Otu012 43.93690 -1.27908 0.438006 -2.92024 0.00349767 0.00983306
res_sig2 = cbind(as(res_sig2, "data.frame"), as(tax_table(microSub)[rownames(res_sig2), ], "matrix"))
ggplot(res_sig2, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
geom_jitter(size=5, width = 0.2) + theme(text = element_text(size=16)) +
theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
# compare 10%Fat and Oleate
res3 = results(ds, contrast=c("condition", "10%Fat", "Oleate"), alpha=alpha)
res3 = res3[order(res3$padj, na.last=NA), ]
res_sig3 = res3[(res3$padj < alpha), ]
res_sig3
## log2 fold change (MLE): condition 10%Fat vs Oleate
## Wald test p-value: condition 10%Fat vs Oleate
## DataFrame with 34 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## Otu007 92.6627 -9.43646 0.542183 -17.40455 7.62021e-68 1.20399e-65
## Otu035 14.4848 -7.15821 0.783378 -9.13761 6.38477e-20 5.04397e-18
## Otu036 13.0374 -6.40371 0.707215 -9.05484 1.36771e-19 7.20329e-18
## Otu033 12.2718 -6.20581 0.720641 -8.61152 7.21001e-18 2.84795e-16
## Otu023 20.7554 -4.48379 0.549579 -8.15858 3.38974e-16 1.07116e-14
## ... ... ... ... ... ... ...
## Otu029 21.614602 -2.79377 0.810938 -3.44510 0.000570838 0.00300641
## Otu030 16.043388 1.23116 0.375613 3.27775 0.001046392 0.00516656
## Otu068 2.640345 -3.05941 0.933212 -3.27836 0.001044109 0.00516656
## Otu125 0.721507 -2.90713 0.906565 -3.20676 0.001342406 0.00642728
## Otu019 21.243604 -2.03839 0.663115 -3.07396 0.002112387 0.00981639
res_sig3 = cbind(as(res_sig3, "data.frame"), as(tax_table(microSubt)[rownames(res_sig3), ], "matrix"))
ggplot(res_sig3, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
geom_jitter(size=5, width = 0.2) + theme(text = element_text(size=16)) +
theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
## Plot heatmap
plot_taxa_heatmap(microSubt,
subset.top = 20,
VariableA = "condition",
heatcolors = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100),
transformation = "log10"
)
## Top 20 OTUs selected
## log10, if zeros in data then log10(1+x) will be used
## First top taxa were selected and
## then abundances tranformed to log10(1+X)
## Warning in transform(phyobj1, "log10"): OTU table contains zeroes. Using log10(1
## + x) transform.
## Venn Diagram
pseq<-microSubt
table(meta(pseq)$condition, useNA = "always")
##
## 10%Fat 60%Fat Oleate <NA>
## 21 14 21 0
# convert to relative abundances
pseq.rel <- microbiome::transform(pseq, "compositional")
disease_states <- unique(as.character(meta(pseq.rel)$condition))
print(disease_states)
## [1] "10%Fat" "Oleate" "60%Fat"
#Write a for loop to go through each of the disease_states one by one and combine identified core taxa into a list.
list_core <- c() # an empty object to store information
for (n in disease_states){ # for each variable n in DiseaseState
#print(paste0("Identifying Core Taxa for ", n))
ps.sub <- subset_samples(pseq.rel, condition == n) # Choose sample from DiseaseState by n
core_m <- core_members(ps.sub, # ps.sub is phyloseq selected with only samples from g
detection = 0.001, # 0.001 in atleast 90% samples
prevalence = 0.75)
print(paste0("No. of core taxa in ", n, " : ", length(core_m))) # print core taxa identified in each DiseaseState.
list_core[[n]] <- core_m # add to a list core taxa for each group.
#print(list_core)
}
## [1] "No. of core taxa in 10%Fat : 38"
## [1] "No. of core taxa in Oleate : 39"
## [1] "No. of core taxa in 60%Fat : 38"
print(list_core)
## $`10%Fat`
## [1] "Otu046" "Otu015" "Otu011" "Otu037" "Otu030" "Otu045" "Otu005" "Otu001"
## [9] "Otu012" "Otu052" "Otu042" "Otu034" "Otu028" "Otu027" "Otu014" "Otu061"
## [17] "Otu002" "Otu075" "Otu013" "Otu017" "Otu016" "Otu025" "Otu003" "Otu004"
## [25] "Otu040" "Otu020" "Otu026" "Otu048" "Otu059" "Otu024" "Otu018" "Otu031"
## [33] "Otu009" "Otu021" "Otu022" "Otu053" "Otu051" "Otu054"
##
## $Oleate
## [1] "Otu046" "Otu015" "Otu011" "Otu037" "Otu074" "Otu045" "Otu005" "Otu001"
## [9] "Otu012" "Otu052" "Otu042" "Otu034" "Otu028" "Otu027" "Otu023" "Otu014"
## [17] "Otu002" "Otu013" "Otu062" "Otu039" "Otu016" "Otu025" "Otu003" "Otu004"
## [25] "Otu008" "Otu044" "Otu041" "Otu040" "Otu020" "Otu026" "Otu006" "Otu059"
## [33] "Otu024" "Otu031" "Otu009" "Otu055" "Otu021" "Otu054" "Otu007"
##
## $`60%Fat`
## [1] "Otu015" "Otu011" "Otu010" "Otu037" "Otu030" "Otu045" "Otu005" "Otu001"
## [9] "Otu012" "Otu033" "Otu029" "Otu034" "Otu023" "Otu014" "Otu076" "Otu002"
## [17] "Otu013" "Otu017" "Otu035" "Otu039" "Otu025" "Otu003" "Otu008" "Otu044"
## [25] "Otu041" "Otu020" "Otu026" "Otu036" "Otu006" "Otu019" "Otu048" "Otu024"
## [33] "Otu018" "Otu078" "Otu009" "Otu021" "Otu050" "Otu007"
# Specify colors and plot venn
# supplying colors in the order they appear in list_core
mycols <- c("brown3", "grey50", "steelblue")
venn.diagram(list_core,fill = mycols, filename = '#14_venn_diagramm.png',
output=TRUE)
## [1] 1
pa <- aggregate_taxa(microSubt, "Genus")
top_four <- top_taxa(pa, 4)
top_four
## [1] "Lachnospiraceae" "Lactococcus" "Muribaculaceae" "Akkermansia"
mycols <- c("brown3", "steelblue", "grey50")
p <- plot_listed_taxa(pa, top_four,
group= "condition",
group.order = c("10%Fat","Oleate","60%Fat"),
group.colors = mycols,
add.violin = TRUE,
violin.opacity = 0.3,
dot.opacity = 0.25,
box.opacity = 0.25,
panel.arrange= "wrap")
## An additonal column Sam_rep with sample names is created for reference purpose
comps <- make_pairs(sample_data(pa)$condition)
p <- p + stat_compare_means(
comparisons = comps,
label = "p.format",
tip.length = 0.05,
method = "wilcox.test")
print(p + ylab("Relative abundance") + scale_y_continuous(labels = scales::percent))
## Get taxa summary by group(s)
p0 <- microSubt
p0.gen <- aggregate_taxa(microSubt,"Genus")
x.d <- dominant_taxa(p0,level = "Genus", group="condition")
head(x.d$dominant_overview)
## # A tibble: 6 x 5
## # Groups: condition [3]
## condition dominant_taxa n rel.freq rel.freq.pct
## <fct> <chr> <int> <dbl> <chr>
## 1 10%Fat Lachnospiraceae 14 66.7 67%
## 2 60%Fat Lachnospiraceae 9 64.3 64%
## 3 Oleate Lachnospiraceae 9 42.9 43%
## 4 10%Fat Lactococcus 6 28.6 29%
## 5 Oleate Lactococcus 6 28.6 29%
## 6 Oleate Akkermansia 4 19 19%
tx.sum1 <- taxa_summary(p0, "Phylum")
## Data provided is not compositional
## will first transform
p0.rel <- microbiome::transform(p0, "compositional")
grp_abund <- get_group_abundances(p0.rel,
level = "Phylum",
group="condition",
transform = "compositional")
## An additonal column Sam_rep with sample names is created for reference purpose
mycols <- c("brown3", "steelblue","grey50")
# clean names
grp_abund$OTUID <- gsub("p__", "",grp_abund$OTUID)
grp_abund$OTUID <- ifelse(grp_abund$OTUID == "",
"Unclassified", grp_abund$OTUID)
mean.plot <- grp_abund %>% # input data
ggplot(aes(x= reorder(OTUID, mean_abundance), # reroder based on mean abundance
y= mean_abundance,
fill=condition)) + # x and y axis
geom_bar( stat = "identity",
position=position_dodge()) +
scale_fill_manual("condition", values=mycols) + # manually specify colors
theme_bw() + # add a widely used ggplot2 theme
ylab("Mean Relative Abundance") + # label y axis
xlab("Phylum") + # label x axis
coord_flip() # rotate plot
mean.plot
## Plotting density of reads per sample
p0 <- microSubt
print_ps(p0)
## 01] ntaxa = 491
## 02] nsamples = 56
## 03] nsamplesvariables = 4
## 04] nranks = 6
## 05] Min. number of reads = 2766
## 06] Max. number of reads = 2766
## 07] Total number of reads = 154896
## 08] Average number of reads = 2766
## 09] Median number of reads = 2766
## 10] Sparsity = 0.792588012801862
## 11] Number of singletons = 179
## 12] % of taxa that are singletons
## (i.e. exactly one read detected across all samples) = 36.4562118126273
kable(head(tax_table(p0)))
| Kingdom | Phylum | Class | Order | Family | Genus | |
|---|---|---|---|---|---|---|
| Otu086 | Bacteria | Actinobacteriota | Actinobacteria | Bifidobacteriales | Bifidobacteriaceae | Bifidobacterium |
| Otu166 | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | uncultured |
| Otu201 | Bacteria | Actinobacteriota | Coriobacteriia | Coriobacteriales | Coriobacteriales_unclassified | Coriobacteriales |
| Otu170 | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | GCA-900066575 |
| Otu125 | Bacteria | Firmicutes | Clostridia | Oscillospirales | Ruminococcaceae | Incertae |
| Otu512 | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospiraceae |
# reduce size for example
ps0 <- core(p0, detection = 10, prevalence = 20 / 100)
# Add a prefix to taxa labels
ps0.f2 <- format_to_besthit(ps0, prefix = "MyBug-")
kable(head(tax_table(ps0.f2))[3:6])
| Domain | Phylum | Class | Order | Family | Genus | best_hit | |
|---|---|---|---|---|---|---|---|
| MyBug-Otu011:Lachnospiraceae | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospiraceae | MyBug-Otu011:Lachnospiraceae |
| MyBug-Otu010:Turicibacter | Bacteria | Firmicutes | Bacilli | Erysipelotrichales | Erysipelotrichaceae | Turicibacter | MyBug-Otu010:Turicibacter |
| MyBug-Otu038:Lachnospiraceae | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospiraceae | MyBug-Otu038:Lachnospiraceae |
| MyBug-Otu037:Oscillibacter | Bacteria | Firmicutes | Clostridia | Oscillospirales | Oscillospiraceae | Oscillibacter | MyBug-Otu037:Oscillibacter |
p <- plot_read_distribution(ps0.f2, groups = "condition",
plot.type = "density")
print(p + theme_biome_utils())
ps0 <- core(p0, detection = 10, prevalence = 20 / 100)
pseq_df <- phy_to_ldf(ps0, transform.counts = NULL)
## An additonal column Sam_rep with sample names is created for reference purpose
kable(head(pseq_df))
| OTUID | Kingdom | Phylum | Class | Order | Family | Genus | Sam_rep | Abundance | group | site | mouse | condition |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Otu046 | Bacteria | Firmicutes | Clostridia | Oscillospirales | Ruminococcaceae | Ruminococcaceae | ce4406 | 34 | ce4406 | cecum | 4406 | 10%Fat |
| Otu015 | Bacteria | Firmicutes | Clostridia | Oscillospirales | Oscillospiraceae | Colidextribacter | ce4406 | 76 | ce4406 | cecum | 4406 | 10%Fat |
| Otu011 | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospiraceae | ce4406 | 132 | ce4406 | cecum | 4406 | 10%Fat |
| Otu010 | Bacteria | Firmicutes | Bacilli | Erysipelotrichales | Erysipelotrichaceae | Turicibacter | ce4406 | 0 | ce4406 | cecum | 4406 | 10%Fat |
| Otu038 | Bacteria | Firmicutes | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospiraceae | ce4406 | 0 | ce4406 | cecum | 4406 | 10%Fat |
| Otu037 | Bacteria | Firmicutes | Clostridia | Oscillospirales | Oscillospiraceae | Oscillibacter | ce4406 | 27 | ce4406 | cecum | 4406 | 10%Fat |
# Plot Density of phyla
pseq <- microSubt
# check 10%Fat
control_ps <- subset_samples(pseq, condition=="10%Fat")
p_hc <- taxa_distribution(control_ps) +
theme_biome_utils() +
labs(title = "10%Fat")
# check Oleate
oleate_ps <- subset_samples(pseq, condition=="Oleate")
p_oleate <- taxa_distribution(oleate_ps) +
theme_biome_utils() +
labs(title = "Oleate")
# check 60%Fat
HFD_ps <- subset_samples(pseq, condition=="60%Fat")
p_hfd <- taxa_distribution(HFD_ps) +
theme_biome_utils() +
labs(title = "60%Fat")
# harnessing the power of patchwork
p_hc / p_oleate / p_hfd + plot_layout(guides = "collect") +
plot_annotation(tag_levels = "A")
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 210 rows containing non-finite values (stat_density).
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 184 rows containing non-finite values (stat_density).
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 203 rows containing non-finite values (stat_density).
## Filtering the dataset
# Select specific taxa
microFirm <- subset_taxa(microSubt, Phylum %in% "Firmicutes")
random_tree2 = rtree(ntaxa(microFirm), rooted=TRUE, tip.label=taxa_names(microFirm))
plot(random_tree2)
microFirm = merge_phyloseq(microFirm, random_tree2)
microFirm
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 426 taxa and 56 samples ]
## sample_data() Sample Data: [ 56 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 426 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 426 tips and 1 internal nodes ]
# Subset samples by site or condition
microFecal <- subset_samples(microSubt, site=="fecal")
microFecal
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 491 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 491 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 491 tips and 490 internal nodes ]
# if you needed to remove candidate outliers, can use the below to remove sample e.g. f4406
microSubt.1 <- prune_samples(sample_names(microSubt) != "f4406", microSubt)
x<-psmelt(microSubt)
head(x)
## OTU Sample Abundance group site mouse condition Kingdom Phylum
## 11 Otu001 f4424 1571 f4424 fecal 4424 60%Fat Bacteria Firmicutes
## 50 Otu001 f4412 1220 f4412 fecal 4412 10%Fat Bacteria Firmicutes
## 24 Otu001 f4418 1196 f4418 fecal 4418 Oleate Bacteria Firmicutes
## 46 Otu001 f4407 1161 f4407 fecal 4407 10%Fat Bacteria Firmicutes
## 9 Otu001 f4406 933 f4406 fecal 4406 10%Fat Bacteria Firmicutes
## 53 Otu001 f4411 929 f4411 fecal 4411 10%Fat Bacteria Firmicutes
## Class Order Family Genus
## 11 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 50 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 24 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 46 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 9 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 53 Bacilli Lactobacillales Streptococcaceae Lactococcus
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 grid stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.1 ggpubr_0.4.0
## [3] RColorBrewer_1.1-2 viridis_0.5.1
## [5] viridisLite_0.3.0 microbiomeutilities_1.00.15
## [7] dplyr_1.0.4 microbiome_1.12.0
## [9] VennDiagram_1.6.20 futile.logger_1.4.3
## [11] microbial_0.0.17 DESeq2_1.30.1
## [13] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [15] MatrixGenerics_1.2.1 matrixStats_0.58.0
## [17] GenomicRanges_1.42.0 GenomeInfoDb_1.26.2
## [19] IRanges_2.24.1 S4Vectors_0.28.1
## [21] BiocGenerics_0.36.0 ape_5.4-1
## [23] ellipse_0.4.2 clustsig_1.1
## [25] knitr_1.31 directlabels_2021.1.13
## [27] vegan_2.5-7 lattice_0.20-41
## [29] permute_0.9-5 plyr_1.8.6
## [31] ggplot2_3.3.3 phyloseq_1.34.0
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 fastmatch_1.1-0
## [4] igraph_1.2.6 splines_4.0.3 BiocParallel_1.24.1
## [7] digest_0.6.27 foreach_1.5.1 htmltools_0.5.1.1
## [10] fansi_0.4.2 magrittr_2.0.1 memoise_2.0.0
## [13] cluster_2.1.1 DECIPHER_2.18.1 openxlsx_4.2.3
## [16] limma_3.46.0 Biostrings_2.58.0 annotate_1.68.0
## [19] RcppParallel_5.0.3 prettyunits_1.1.1 jpeg_0.1-8.1
## [22] colorspace_2.0-0 ggrepel_0.9.1 blob_1.2.1
## [25] haven_2.3.1 xfun_0.21 crayon_1.4.1
## [28] RCurl_1.98-1.2 jsonlite_1.7.2 genefilter_1.72.1
## [31] dada2_1.18.0 phangorn_2.5.5 survival_3.2-7
## [34] iterators_1.0.13 glue_1.4.2 gtable_0.3.0
## [37] zlibbioc_1.36.0 XVector_0.30.0 DelayedArray_0.16.1
## [40] car_3.0-10 Rhdf5lib_1.12.1 abind_1.4-5
## [43] scales_1.1.1 pheatmap_1.0.12 futile.options_1.0.1
## [46] DBI_1.1.1 edgeR_3.32.1 rstatix_0.7.0
## [49] Rcpp_1.0.6 xtable_1.8-4 progress_1.2.2
## [52] foreign_0.8-81 bit_4.0.4 httr_1.4.2
## [55] ellipsis_0.3.1 farver_2.0.3 pkgconfig_2.0.3
## [58] XML_3.99-0.5 sass_0.3.1 locfit_1.5-9.4
## [61] utf8_1.1.4 labeling_0.4.2 tidyselect_1.1.0
## [64] rlang_0.4.10 reshape2_1.4.4 AnnotationDbi_1.52.0
## [67] munsell_0.5.0 cellranger_1.1.0 tools_4.0.3
## [70] cachem_1.0.4 cli_2.3.1 generics_0.1.0
## [73] RSQLite_2.2.3 ade4_1.7-16 broom_0.7.5
## [76] evaluate_0.14 biomformat_1.18.0 stringr_1.4.0
## [79] fastmap_1.1.0 yaml_2.2.1 bit64_4.0.5
## [82] zip_2.1.1 randomForest_4.6-14 purrr_0.3.4
## [85] nlme_3.1-152 formatR_1.7 rstudioapi_0.13
## [88] compiler_4.0.3 curl_4.3 png_0.1-7
## [91] ggsignif_0.6.1 tibble_3.0.6 geneplotter_1.68.0
## [94] bslib_0.2.4 stringi_1.5.3 highr_0.8
## [97] forcats_0.5.1 Matrix_1.3-2 multtest_2.46.0
## [100] vctrs_0.3.6 pillar_1.5.0 lifecycle_1.0.0
## [103] rhdf5filters_1.2.0 jquerylib_0.1.3 data.table_1.14.0
## [106] bitops_1.0-6 R6_2.5.0 latticeExtra_0.6-29
## [109] hwriter_1.3.2 ShortRead_1.48.0 gridExtra_2.3
## [112] rio_0.5.16 codetools_0.2-18 lambda.r_1.2.4
## [115] MASS_7.3-53.1 assertthat_0.2.1 rhdf5_2.34.0
## [118] withr_2.4.1 GenomicAlignments_1.26.0 Rsamtools_2.6.0
## [121] GenomeInfoDbData_1.2.4 mgcv_1.8-34 hms_1.0.0
## [124] gghalves_0.1.1 quadprog_1.5-8 tidyr_1.1.2
## [127] rmarkdown_2.7 carData_3.0-4 Rtsne_0.15